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ABSTRACT 

A method of adapting smoothed particle hydrodynamics (SPH) with periodic bound- 
ary conditions for use with the special purpose device Grape is presented. Grape 
(GRAvity PipE) solves the Poisson and force equations for an N-body system by di- 
rect summation on a specially designed chip and in addition returns the neighbour list 
for each particle. Due to its design, Grape cannot treat periodic particle distributions 
directly. This limitation of GrapeSPH can be overcome by computing a correction 
force for each particle due to periodicity (Ewald correction) on the host computer 
using a PM-like method. 

This scheme is applied to study the fragmentation process in giant molecular 
clouds. Assuming a pure isothermal model, we follow the dynamical evolution in the 
interior of a molecular cloud starting from an Gaussian initial density distribution to 
the formation of selfgravitating clumps until most of the gas is consumed in these dense 
cores. Despite its simplicity, this model can reproduce some fundamental properties of 
observed molecular clouds, like a clump mass distribution of the form dN/dm oc to", 
with n ~ —1.5. 

Key words: Hydrodynamics - Methods: numerical - ISM: clouds ~ Stars: formation 



1 INTRODUCTION 



In the early 1990's, the first Grape boards became avail- 
able to the scientific community. Grape is an acronym for 
"GRAvity PipE" and is a hardware device to solve the Pois- 
son and force equations for an N-body-system by direct sum- 
mation on a specially designed chip, thus leading to consider- 
able speed-up in computing this type of problems (Sugimoto 
et al. 1990, Ebisuzaki et al. 1993). Besides the gravitational 
forces, Grape also returns the list of neighbours for each 
particle, which makes it attractive for use in combination 
with smoothed particle hydrodynamics, SPH, (for a gen- 
eral review on SPH see Benz 1990 or Monaghan 1992, for 
combination with the hierarchical tree method, see Hern- 
quist & Katz 1989, and for implementing with Grape, see 
Umemura et al. 1993, and Steinmetz 1996). In the last five 
years. Grape boards have been used for studying a variety 
of astrophysical problems, ranging from following the evo- 
lution of binary black holes (Makino et al. 1993a) up to the 
formation of galaxy clusters in a cosmological context (Huss, 
Jain & Steinmetz 1997). A complete overview is given in the 
review article by Makino & Hut (1997). 

Due to its restricted force law. Grape cannot treat pe- 
riodic particle distributions. This limits its applicability to 
strongly self gravitating systems. However, many astrophys- 
ical problems require a correct treatment of periodic bound- 



ary conditions. Simulating the dynamical evolution and frag- 
mentation in the interior of a molecular cloud, a process that 
finally leads to the formation of new stars, requires periodic 
boundaries to prevent the whole object from collapsing to 
the center. The situation is quite similar in cosmological 
large scale structure simulations, where periodic boundaries 
mimic the homogeneity and isotropy of the initial matter dis- 
tribution. In the following, we present a method to overcome 
this limitation and to incorporate fully periodic boundary 
conditions into an N-body or SPH algorithm using Grape. 
The basic idea is to compute a periodic correction force for 
each particle on the host computer, applying a particle-mesh 
(PM) like scheme: We first compute the forces in the isolated 
system using direct summation on Grape, then we assign 
the particle distribution to a mesh and compute the cor- 
rection force for each grid point, by convolution with the 
adequate Green's function in Fourier space. Finally we add 
this correction to each particle in the simulation. The cor- 
rective Green's function can be constructed as the offset be- 
tween the periodic solution (calculated via the Ewald (1921) 
approximation) and the isolated solution on the grid. 

Our approach is thus related to the method described in 
Huss et al. (1997). However, they apply the combined force 
calculation via a periodic PM scheme and direct summation 
on Grape only to the particle distribution in the center of 
their simulation volume. In this region, they compute the 
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dynamical evolution of a galaxy cluster with high resolu- 
tion, whereas the remaining part of the volume is treated 
by a pure PM scheme and provides the tidal torques for the 
cluster in the center. For this reason, they do not have to 
handle possible force misalignment in the border region of 
the simulation cube. In our situation, studying the evolution 
and fragmentation in the interior of molecular clouds, we 
need high resolution and periodicity throughout the whole 
simulation volume. Ensuring numerical stability throughout 
the complete volume requires additional care in combining 
both methods, the PM scheme and Grape. 

This paper is organised as follows: The next section 
gives a brief overview of the main properties of the spe- 
cial purpose hardware Grape; then follows a description 
of the Ewald method which leads to the construction of a 
Green's function for a periodic correction. This Green's func- 
tion is applied in a PM like scheme, to obtain a (corrective) 
force term which is added to the particle forces in Grape- 
SPH. Section ^ investigates the performance of the method 
in selected test cases and suggests ways of optimising the 
algorithm. The next section applies the above to a prelim- 
inary study of the fragmentation process in the interior of 
molecular clouds and the section H finally concludes with a 
summary of this paper. 



2 GRAPE SPECIFICATIONS 

Grape is a special purpose hardware device, which calcu- 
lates the forces and the potential in the gravitational N-body 
problem by direct summation on a specifically designed chip 
with high efficiency (Sugimoto et al. 1990, Ebisuzaki et al. 
1993). We use the currently distributed version, Grape- 
3AF, which contains 8 chips on one board and therefore 
can compute the forces on 8 particles in parallel. The board 
is connected via a standard VME interface to the host com- 
puter, in our case a SUN workstation. C and FORTRAN li- 
braries provide the software interface between the user's pro- 
gram and the board. The computational speed of Grape- 
3AF is approximately 5 Gflops. 

The force law is hardwired to be a Plummer law, 

^ (Iri - rjf + e2)'^/2 

Here i is the index of the particle for which the force is calcu- 
lated and j enumerates the particles which exert the force; 
ti is the gravitational smoothing length of particle i; G and 
rrii and rrij are Newton's constant and the particle masses, 
respectively. To perform the force calculation for a single 
particle, Grape-3AF needs A*' -f- 19 clock cycles. To achieve 
a high speed, concessions in the accuracy of the force calcu- 
lations had to be made: Grape internally works with a 20 
bit fixed point number format for particle positions, with a 
56 bit fixed point number format for the forces and a 14 bit 
logarithmic number format for the masses (Okumura et al. 
1993). Conversion to and from that internal number repre- 
sentation is handled by the interface software. The number 
format limits the spacial resolution in a simulation and con- 
strains the force accuracy. However, for coUisionless N-body 
systems, the forces on a single particle typically need not 
be known better than up to an error of about one percent. 



In that respect, Grape is comparable to the widely used 
Treecode schemes (e.g. Barnes & Hut 1986). Besides the 
low precision Grape board, there is a high precision version 
available. Harp, which was specifically designed for treat- 
ment of collision dominated systems (Makino, Kokubo and 
Taiji 1993). 

In addition to integrate the time evolution of an N-body 
system by direct summation. Grape has been incorporated 
into a Treecode algorithm (Makino & Funato 1993) and, 
more recently, into a P3M code (Brieu, Summers & Ostriker 
1995). Besides forces and potential. Grape returns the list 
of neighbours for each particle i within a specified radius hi. 
This feature makes it attractive for use in smoothed particle 
hydrodynamics (Unemura 1993), where gas properties are 
derived by averaging over a localised ensemble of discrete 
particles. For particles near the surface of the integration 
volume, "ghost" particles are created to correctly extend 
the neighbour search beyond the cell border; no forces are 
computed for these particles. For periodic boundaries, the 
ghosts are replicas of particles from the opposite side of the 
cell. Following the strategy described in Steinmetz (1996), 
we have modified our SPH code (Benz 1990, Bate et al. 1995) 
to incorporate Grape and furthermore added the ability to 
handle periodic boundary conditions. 

Note, that there also has been an attempt to treat pe- 
riodic boundaries on a special purpose hardware device. 
Wine, using the Ewald method hardwired on a special 
chip (Fukushige et al. 1993). This work has been success- 
ful. Fukushige et al. (1996) report the development of the 
MD-Grape, which can handle an arbitrary central force 
law, including the Ewald method. 



3 IMPLEMENTING PERIODIC BOUNDARIES 
INTO GRAPE 

The general concept of incorporating periodic boundaries 
into Grape is the following: We use the Ewald (1921) 
method to optain a Green's function that accounts for the 
correction to potential and force for a system transformed 
from having vacuum boundary conditions to periodic ones. 
Grape computes the forces and potential for the isolated 
system. The Green's function is then used in a PM like 
scheme to assign a correction term due to periodicity to 
each particle in the simulation. 

3.1 Periodic Force Correction — The Ewald 
Method 

In 1921, P. Ewald suggested a method to compute the forces 
in a periodic particle distribution (he aimed explicitely to 
compute the potential in atomic lattices in solids). For 
a more recent discussion see Hernquist, Bouchet & Suto 
(1991); for applications to large scale structure simulations 
see Katz, Weinberg & Hernquist (1996), or Dave, Dubinski 
& Hernquist (1997). 

Poisson's equation for a system of N particles which 
are infinitely replicated in all directions with period L is (n 
being an integer vector) 

JV 

V^0(r) = AnGpir) = 4^G ^ ^ m,S{r - r, - nL) (2) 

n j = l 
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and can be solved in general using the appropriate Green's 
function Q: 
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(3.a) 
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which is for the gravitational potential Q{r) = 1/r, or 
^(k) = ATv/k^ in Fourier space. The sum in Eqn. 3 converges 
very slowly, which strongly limits its numerical applicabil- 
ity. However, Ewald (1921) realised, that convergence can 
be improved considerably by splitting the Green's function 
into a short range Qs and a long range part Ql and solving 
the first in real space and the latter one in Fourier space: 



1 - 47r 

Qs{r) = -erfc(ar) ; C7s(k) = -— 
r 



exp 



4a2 



Ql{v) — ierf(Qfr 



e.(k) = -exp - — 



(4.a) 
(4.b) 



Here a is a scaling factor in units of inverse length and erf (x) 
is the error function with erfc(a;) its complement. 



erf (a;) = 



(5) 



With the Green's function split into two parts, Q — Gs + Gl, 
the potential reads as <j){r) = 0s (r) + ^^(r): 
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(6) 



Finally, the force (say exerted onto particle i) is 



F, = F(rO = -m^\/(t>{r,) = Gm, ) m,f(r, -r,), (7) 



with 
f(r) 
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nL 
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erfc (a|r — nL|) + 



H — -=|r — nLj exp(— a |r — nL]' 



(8) 



To achieve good accuracy with reasonable computational 
effort, typical values are a — 2/L, |r — nL| < 3.61/ and 
k = 27rh/L with h being an integer vector with jhp < 10 
(see e.g. Hernquist et al. 1991; note, the second component 
in their Eqn. 2.14b is missing a factor |r — nL|; it is identical 
to our Eqn. ^). 

For a particle pair with separation r, f(r) includes the 
contributions from all L-periodic pairs with identical sepa- 
ration. To obtain the pure periodic correction force fcor(r), 
one has to subtract the direct interaction of the central pair, 
which leads to 



fco..(r) = f(r) 



(9) 



To obtain the periodic correction for the potential of the par- 
ticle pair, (^cor(r), one proceeds similar and again subtracts 
the isolated solution from the Ewald solution (Eqn. ^), 



r(r) = Hr) + 



(10) 



We compute the correction terms fcor and (j)coi for particle 
pairs on a Cartesian grid covering our whole simulation box, 
placing particle 1 in the central node and particle 2 on dif- 
ferent grid points, and obtain so a table of pairwise force 
values. 



3.2 The PM method 

Particle-mesh methods assign particle properties to mesh 
points, solve the interaction equations on the grid and inter- 
polate the solution back onto the particles (see e.g. Hock- 
ney & Eastwood 1988). Typically, the density at each grid 
point is determined from the particle positions and masses 
using CIC ("cloud in cell") or TSC ("triangular shaped 
cloud") schemes, whose assignment functions are triangles 
or quadratic splines, respectively (again Hockney & East- 
wood 1988, chap. 5). For the gravitational N-body prob- 
lem one solves Poisson's equation. Usually this is done in 
Fourier space, since the differential operation there acts as 
a simple multiplication. Therefore, the density distribution 
on the grid is transformed into fc-space using EFT and con- 
volved with the appropriate Green's function to solve for 
the potential. To obtain the forces, we convolve the Fourier 
transform of the density with the Green's function for the 
force^; separately for the x, y and z-component. Finally, in- 
verse FFT returns potential and force at each grid point. 
The alternative for computing the forces in Fourier space 
is, to calculate them in real space as gradients of the mesh- 
defined potential. This would induce additional errors and 
therefore we prefer the first method. 



3.3 Combining the PM Method with Direct 
Summation on GRAPE 

In our application of the PM method, we compute the (pe- 
riodic) correction to the force and the potential and add 
this to the solution for the isolated system, which is ob- 
tained via direct summation on Grape. This procedure en- 
sures proper treatment of periodic boundary conditions. The 
Green's function for the PM scheme can be constructed di- 
rectly in Fourier space (Hockney & Eastwood 1988). How- 
ever, for the correction force and potential, this is rather 
complicated and it is more intuitive and straightforward to 
proceed the following way: We obtain the Green's function 
for each force component as the Fourier transform of the x, y 
and z-component, respectively, of the mesh-defined pairwise 
periodic correction force fcor, defined in Eqn. ^. And for the 
potential correction, the appropriate Green's function is the 
Fourier transform of cjf>cor in Eqn. nol The corrective Green's 



* Strictly speaking, a Green's function is the solution of an equa- 
tion of the type, V^G(i') = — 47n5(r) and its boundary value prob- 
lem. Green's functions are used to solve Poisson's equation for 
the potential. Being the gradient of Poisson's equation, the force 
equation can be solved by the same ansatz. 
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Figure 1. Matrix of pairwise forces - x-component F^; cut through the yz-plane at x=0: a) The periodic contribution (replicated onto 
the whole grid), b) The isolated solution, c) The difference between both, i.e. the correction term. 



function is thus the offset between the Green's function for 
the periodic system C/per and the isolated one Giso, 



Gcoi — ^pcr Gis 



(11) 



Or with other words, fcor and ^cor are the Green's functions 
with the right properties and only have to be transformed 
into Fourier space for convolution with the density. Using 
the Fourier transform of the Ewald forces as periodic Green's 
function was proposed by A. Huss (private communication) 
and it is straightforward to extend this for handling the force 
correction by subtracting the isolated solution. 

Note, in order to obtain the isolated solution in Fourier 
space, one has to double the linear dimensions of the grid 
and zero pad the additional grid points to avoid contami- 
nation from implicitly assumed periodicity (see Press et al. 
1989). For example, to solve Poisson's equation for a cubic 
density field with volume L^, we have to use a grid of size 
{2L)^, assign the density field to one octant of the large box 
and fill the remaining grid points with zero. The Green's 
function, however, has to be defined on the complete grid, 
i.e. G{r) = l/|r| for the potential, with —L < x,y, z < L. On 
the other hand, the periodic Green's function is defined on 
the original grid, with —L/2 < x,y, z < L/2. For alignment 
with the isolated solution, we have to extend the periodic 
solution into the larger cube. To obtain the right correction 
Green's function, e.g. for the force, we replicate the table of 
pairwise Ewald forces (defined on the small cube) into all oc- 
tants of the large cube. Then we subtract the isolated force 
field (defined on the large cube) and transform into Fourier 
space. This procedure is illustrated in Fig. |l[ it plots in 
the yz-plane at a; = 0: a) is the Ewald periodic force com- 
puted on a 32-grid and replicated four time into a 64-grid, 
b) is the isolated solution defined on the large grid and c) is 
the difference, the final correction force. The Green's func- 
tion finally is the Fourier transform of this force matrix. 
Convolution with the zero padded density field results in 
the right force correction. The procedure is analogue for the 
potential. In practice, it is sufficient to compute the Green's 
function once at the beginning of a simulation run and store 
it as a table. 

The total force acting on an individual particle during 
one timestep then stems from the particle system inside the 
simulation cube, computed by direct summation on Grape, 
plus the contribution from an infinite set of periodically mir- 
rored systems, computed via the method described above. 
For smoothed particle hydrodynamics, one also has to add 
pressure and viscous forces. 



We compute the periodic and the isolated solution on 
a grid, subtract the latter from the first and add the iso- 
lated forces calculated with Grape. The two isolated solu- 
tions cancel out and one can ask the question, what have we 
gained? For a system with more or less homogeneous den- 
sity a pure PM scheme is sufficient. However, the advantage 
of using Grape is evident when computing strongly struc- 
tured systems. Unlike in PM schemes, the spacial resolution 
with Grape is not limited to a given cell size, but adapts to 
the density distribution due to its Lagrangian nature. It is 
limited only by the gravitational smoothing length, or equiv- 
alently, by the choice of the minimum timestep. Since the 
potential of strong density peaks is dominated by self grav- 
ity and the influence of periodic boundaries (and thus of 
the Ewald correction) becomes weak, it is sufficient to com- 
pute this correction term on a relatively coarse grid which 
keeps the additional computational cost low. For a relatively 
smooth density distribution, a relatively wide mesh is suffi- 
cient anyhow. The scheme proposed here unites both, high 
resolution with Grape and the periodicity of a PM scheme. 
In addition, applying a Courant-Friedrichs-Lewy like crite- 
rion, we typically do not have to compute the FFT at each 
smallest timestep, but can use stored correction values from 
the previous call. This reduces the computational expense 
further. 



4 OPTIMIZATION OF THE SCHEME AND 
PERFORMANCE IN TEST CASES 

The contributions to the total acceleration of one particle are 
computed by two distinct methods: by direct summation on 
the Grape board for the isolated solution and by applying a 
particle-mesh scheme for the periodic correction term. Com- 
pared to the host computer, the Grape c hips hereby have 
only limited numerical accuracy (see Sect. 4.1). The spacial 
resolution of the PM scheme is limited by the number of gri d 
zones and the choice of the assignment function (Sect. 4.2). 
All this may lead to spurious residuals when combining the 
forces; some terms might not cancel out completely. This 
"misalignment" can be minimised by randomly shifting the 
center of the simulation cube through the periodic particle 
distribution. Then on average the force residuals cancel out 
and numerical stability is increased (Sect. 1.3). 
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Figure 2. Influence of the gravitational smoothing length e. For 
all e < 0.001 the force errors are equally distributed and minimal 
(thick dashed line). The errors grow with increasing e (e = 0.005 
— solid line, e = 0.01 dashed line). These values are computed 
for a distribution of 32'^ particles on a regular grid with periodic 
boundaries. The force correction is computed on a 64^-mesh. The 
particles are centered on grid points to be independent of the 
assignment scheme. 

4.1 Gravitational Smoothing Length and 
Numerical Accuracy of GRAPE 

Since Grape is used to calculate the force contribution from 
the isolated system, we also call the board at the start of the 
simulation, when setting up the force correction table. We 
subtract the isolated solution from the periodic Ewald solu- 
tion for each particle pair on a grid. For our test problems, 
we use a fixed gravitational smoothing length e throughout 
the whole computation and also for calculating the periodic 
Ewald forces. As always in coUisionless N-body simulations, 
the choice of e is a trade-off between minimising 2-body re- 
laxation (large e) and spacial resolution (small e). For the 
simulations presented here, we typically adopt values be- 
tween e = 0.01 and 0.001, with the total size of the simula- 
tion cube being [—1,-1-1]^. The infiuence of e is demonstrated 
in Fig. ^; it plots the force in x-direction. Fx, in dimension- 
less units for 32'' particles which are distributed on a regular 
grid. Assuming periodic boundary conditions, this distribu- 
tion is - in principle - force free. For all values of e < 0.001 
the distribution of force errors is equal and determined by 
the numerical accuracy of the Grape chips. For larger e, 
the force contributions from the two methods get increas- 
ingly misaligned. 

4.2 Grid Resolution and Assignment Function 

Another important factor is the resolution of the grid. The 
more grid zones are used, the smaller the wavelengths that 
can be resolved by the PM scheme. However, the number 
of CPU cycles increases linearly with the number of grid 
zones. Again one has to find an compromise between accu- 
racy and computational speed, depending on the problem to 
be solved. The errors furthermore depend on the adopted as- 
signment function. To obtain the force correction, one has to 
interpolate the particle distribution onto a grid, solve Pois- 
son's equation for the density field to obtain forces and as- 
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Figure 3. Comparison of the relative force errors dF/F for dif- 
ferent grid resolutions: 16'^ (dotted line), 32'' (dashed line) and 
64"' (continuous line). As reference, the error distribution stem- 
ming solely from Grape is plotted as thick dashed line (analogue 
to Fig. |). 

sign these forces back onto the particles. One possibility is 
to assign all particles to the cell they are located in, the 
nearest grid point scheme, NGP. Another way is CIC (cloud 
in cell), which uses a boxlike cloud shape to distribute each 
particle into eight neighbouring cells (in three dimensions). 
TSC (triangular shaped cloud) distributes the particle mass 
into 27 cells (for more details consult Hockney & Eastwood 
1988). For the problems we study here, we consider NGP as 
too coarse, higher order schemes on the other hand smear 
out the particle distribution and thus limit the spatial reso- 
lution. We adopt the CIC scheme in our calculations. Using 
CIC, we plot in Fig. ^the errors in F^ for the distribution 
of 32"' particles on a regular periodic grid, but now shifted 
by dx = 0.015625, i.e. half of the ceU size of a 64^-grid. As 
expected for the finer meshes, the interpolation is more accu- 
rate and the force errors are smaller. As a reference, we have 
plotted the intrinsic error contribution coming from Grape 
alone (the thick long dashed line). Note that the errors are 
not evenly distributed amongst the particles; the errors are 
small and close to zero in the interior of the simulated cube, 
whereas particles at the border have errors by one or two 
orders of magnitudes larger (dotted and dashed lines in the 
range —2 < logj,, dF/F < —1). Note as well, that each grid 
is in fact a factor of 2"' larger to account for the zero padding 
necessary for the isolated solution, i.e. the 64^-grid is in fact 
a 128^-grid. 

4.3 The Influence of Random Shifts 

We can increase the stability of the code by minimising the 
influence of the force errors for the border cells. One way to 
do this is to ensure, that these border cells do not always 
contain the same group of particles: We apply a random shift 
of the whole simulation area. We can think of our simulation 
box as a window onto an infinite periodic particle distribu- 
tion. Since this distribution is infinite and periodic we are 
free to choose any region in it, as long as it contains the 
periodicity of the whole distribution. By randomly shifting 
the center of the simulation box we prevent the errors made 
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Figure 4. Average particle displacement as function of time for 
the system integrated with (squares) and without (diamonds) ran- 
dom shift. 



in the border cells to add up coherently for the same par- 
ticles. On average they cancel out. The stabilising effect of 
this method is demonstrated in Fig. ^ It compares the time 
evolution of the average particle displacement (due to force 
errors) for the standard particle distribution (particles on a 
regular 32'^ mesh) for the two cases with and without apply- 
ing random shift. Clearly the system computed without shift 
degrades much faster than the other one. In fact we plot the 
pathological case, where the particles are located exactly be- 
tween two mesh points and therefore the assignment error 
is maximal. This system collapses due to this errors on a 
timescale of about 5 free-fall timesQ For a system where the 
particles are placed exactly on grid points the assignment 
errors are minimal. Such a system is more stable without 
random shift. However, in a physically meaningful context 
the particles are far from being placed exactly on grid points 
and random shift is an important tool of stabilising the code. 




0.10 
B 



Figure 5. Power spectrum P{k) and 2-point correlation func- 
tion ^{R) for the collapse of a random fluctuation flcld in pe- 
riodic boundaries. The solid lines describe the calculation done 
with Grape and the dashed lines with a Treecode scheme. Both 
systems evolve almost identically. Only for small to intermediate 
particle separations, the difference in the force calculation in both 
schemes leads to a deviation in the correlation and in the power 
spectrum. 



4.4 Comparison with pure Ewald Method 

A further test of the performance of this new method is, 
to compare it with a Treecode scheme. In this case, the 
Ewald method can be directly implemented into the force 
summation, as described by Hernquist et al. (1991). To do 
so, we distribute 32"^ particles randomly within the simula- 
tion cube to get a homogeneous distribution with Poisson 
fluctuations. We chose the initial kinetic energy to be zero 
and follow the fragmentation and collapse of the system. 
As expected for this purely gravitational N-body system, 
density peaks start to grow to form a network of filaments 
and knots, similar to those known from cosmological N-body 
simulations - note however, that we do not include a Hubble 
expansion term. These knots finally merge to form one big 
clump. The evolution of the power spectrum P{k) and of 

t The regular infinite system does not collapse; the free-fall time 
is per default infinite. Here and in the following free-fall time 
means the interval, an isolated simulation cube would need to 
collapse. 



the correlation function (,{R) is illustrated in Fig. q. Con- 
sidering the differences in the force calculation method, the 
evolution of the two systems agree remarkably well. 



5 FRAGMENTATION OF MOLECULAR 
CLOUDS 

As a first application of the method presented above we de- 
scribe the simulation of the dynamical evolution and frag- 
mentation processes in the interior of giant molecular clouds. 

Molecular clouds are highly structured: Observations 
reveal a hierarchy of filaments and sheets, clumps and sub- 
clumps, ranging from scales of the size of the whole cloud 
(see e.g. Bally 1996) down to the smallest objects resolved 
by todays telescopes (Wiseman & Ho 1996). Density-size or 
linewidth-size relations indicate that typically a cloud is to 
a large extent supported by "turbulent motion"; the mea- 
sured linewidths exceed the thermal line broadening by far. 
Observations of Zeeman splitting and polarisation indicate 
the presence of magnetic fields. Furthermore feedback mech- 
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Figure 6. Evolution of a Gaussian density field with P(k) oc 1/fc^: Initial condition (left side) and after one free-fall time (to the right). 
The density is normalised equally and scaled logarithmically in both images. The simulation was done with fully periodic GrapeSPH 
code using 100 000 particles. 



anisms from newly formed stars, outflows, stellar winds, ion- 
ization fronts and finally supernovae, produce shells and 
bubbles and deposit huge amounts of energy and momen- 
tum into the interstellar medium. The physical processes 
in molecular clouds are extremely diverse. A comprehensive 
overview can be found in "Protostars and Planets III" , eds. 
Levy & Lunine, 1993) 

However, this is the environment in which stars form. 
To assess the problem of fragmentation and clump forma- 
tion that leads to star formation, we start in the most basic 
way: We follow the dynamical evolution of isothermal gas 
in a region in the interior of a molecular cloud, starting 
from an initial density distribution to the formation of self- 
gravitating clumps. We solve the hydrodynamical equations 
using the SPH with Grape and implement periodic bound- 
ary conditions as described above to prevent global collapse. 
Ignoring additional physical effects, we study the interplay 
between gravity and gas pressure, which by itself will pro- 
duce hierarchical filamentary structure, as was derived by 
de Vega, Sanchez & Combes (1996). Figure ^ describes one 
simulation: The left image is the projection of the initial 
density field, that was assumed to be Gaussian with a 1/fc^ 
power spectrum. The right image depicts the same field after 
one free-fall time^ The initial Gaussian fluctuations evolved 
into a system of filaments and knots, some of which contain 
collapsing cores. The initial conditions in Fig. ^ where cho- 
sen to be very cold; the total mass in the system exceeds the 
Jeans mass by far and a large number of initial density peaks 
are selfgravitating and begin to collapse. Gravity dominates 
strongly over pressure forces. 

Figure ^ shows the mass distribution of the identified 



t Again, free-fall time means the time interval an isolated system 
would need to collapse. 



clumps. We divide the mass range into equal logarithmic 
intervals and plot the number of clumps identified, divided 
by the interval. For masses below m ~ 0.001, our clump 
finding algorithm deteriorates and the sampling gets incom- 
plete. With the chosen initial conditions and parameters, the 
mass distribution follows a power law, dN/dm oc m", with 
n ~ —1.5, which is indicated by the dashed line. This is in 
agreement with the values observed in molecular clouds (see 
e.g. Blitz 1993) although the observations exhibit a huge 
scatter. These results of our study of the fragmentation pro- 
cess in molecular clouds are still preliminary; more details 
will be discussed in a subsequent paper (Klessen & Burkert 
1997). 

Certainly this treatment is very coarse. It neglects the 
presence of magnetic fields and the input of energy and mo- 
mentum by young stars. Future numerical studies of the evo- 
lution of interstellar medium and the formation of stars have 
to take these processes into account. However, the simple 
isothermal model presented here is able to reproduce some 
of the observed properties of giant molecular clouds and star 
forming regions. 



6 SUMMARY 

We have presented a method of how to incorporate peri- 
odic boundary conditions into numerical simulations using 
the special hardware device Grape. We compute a periodic 
correction force onto each particle applying a particle-mesh 
like scheme: The particle forces in the isolated system are 
calculated via direct summation on the Grape board. Then 
we assign the particle distribution to a mesh and compute 
the periodic correction force for each grid point by convolu- 
tion with the appropriate Green's function in Fourier space. 
These forces are assigned back onto the individual parti- 
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Figure 7. Mass distribution of the identified clumps. The dashed 
line indicates the slope dN/dm oc m~^'^. The total mass in the 
system is m = 1; for clump masses m < O.OOf the sample gets 
incomplete. 



cles. This method unites the advantages of Grape as being 
Lagrangian and offering good spacial resolution necessary 
for treating highly structured density distributions with the 
advantage of particle-mesh schemes as being intrinsically pe- 
riodic. The potential of strong density peaks is dominated 
by the local mass accumulation. In such a region the influ- 
ence of the global properties of the system is week and it is 
therefore sufficient to compute the correction term on rela- 
tively coarse grid. This keeps the additional computational 
cost low. 

The method proposed here has proved to be numerically 
stable and inexpensive in terms of computer time. Thus the 
special purpose hardware device Grape can be applied suc- 
cessfully to astrophysical problems that require the correct 
treatment of periodic boundary conditions. As an example 
we study the dynamical evolution and fragmentation process 
in the interior of molecular clouds. 
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